function f = covar(theta2,resid,sj,v,x2,p,QI,ind_market,A1,W,d_d,ztilde,xtilde,deltaold,sales_tax)

% this function computes the var-covar matrix; 

N=length(resid)/2;
k2=length(theta2); % number of nonlinear parameters
z_D=ztilde(1:N,1:d_d); % demand instruments 
z_S=ztilde(N+1:end,d_d+1:end); % supply instruments 
d_S=size(z_S,2);
xi=resid(1:N,1);
omega=resid(N+1:end,1);


h=1e-8;
H=eye(k2,k2)*h;
[ai, bi]=computey(theta2+H(1,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
[af, bf]=computey(theta2-H(1,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax); 
K=([ai;bi]-[af;bf])/(2*h); 
if k2>1;
    for i=2:k2;
        [Ai,Bi]=computey(theta2+H(i,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
        [Af,Bf]=computey(theta2-H(i,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
        K=[K ([Ai;Bi]-[Af;Bf])/(2*h)];
    end;
end;

Gamma=(ztilde'*[-xtilde K]);

f=inv(Gamma'*W*Gamma); 